Visualizing Spatial Data in R

Kristin Broms, Neptune & Co.

October 18, 2017

Outline

As you can see, “spatial data” has come to mean much more than creating a map! Visualizing spatial data is an active area of software development, with many exciting packages and functions becoming available on a regular basis.

Static map- shapefiles

library(rgdal)
# streams <- readOGR(dsn = "../data/R_GIS_Layers/", 
#                    layer = "Cove_Drainage_WGS84")
summary(streams)
## Object of class SpatialLinesDataFrame
## Coordinates:
##          min        max
## x -109.35383 -109.12624
## y   36.46573   36.94164
## Is projected: FALSE 
## proj4string :
## [+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
## Data attributes:
##     OBJECTID          ComID               GNIS_ID       GNIS_Name 
##  Min.   :     0   Min.   :        0   00003435: 2   Cove Wash: 2  
##  1st Qu.: 32992   1st Qu.:103664065   NA's    :47   NA's     :47  
##  Median : 65593   Median :103665843                               
##  Mean   : 68116   Mean   :101439114                               
##  3rd Qu.:105309   3rd Qu.:126660550                               
##  Max.   :136421   Max.   :126669825                               
##                                                                   
##           ReachCode      FType         FCode        Source  
##  14080204001554: 2   Min.   :460   Min.   :    0   DOQQ: 1  
##  14080204001982: 2   1st Qu.:460   1st Qu.:46003   NHDH:48  
##  14080105001404: 1   Median :460   Median :46003            
##  14080105001406: 1   Mean   :462   Mean   :45264            
##  14080105001455: 1   3rd Qu.:460   3rd Qu.:46003            
##  (Other)       :41   Max.   :558   Max.   :55800            
##  NA's          : 1                                          
##                Name   
##  Middle 1        : 2  
##  Cove Wash       : 1  
##  Cove Wash Middle: 1  
##  Cove Wash North : 1  
##  Cove Wash South : 1  
##  (Other)         :17  
##  NA's            :26

Static map- shapefiles

plot(streams, col="blue")

## note that plots using base R may be distorted.

Static map- shapefiles

library(ggplot2)
streams_gg <- fortify(streams)  ## or use broom::tidy()
# str(streams_gg)
ggplot() + 
  geom_path(data = streams_gg, 
            aes(x = long, y = lat, group = group), 
                color = "blue", size = 1.5)

## note that without a coordinates correction, the plot may be distorted.

Static map- satellite maps

library(ggmap)
# Might need older version of ggplot2:
# devtools::install_github("hadley/ggplot2@v2.2.0")
# library(ggplot2)
myLocation <- c(lon = -109.2279, lat = 36.545) 
myMap <- get_map(location = myLocation, source = "google", 
                 maptype = "satellite", crop = FALSE, zoom = 13)
g <- ggmap(myMap, darken = c(0.25, "white"))
g + 
  geom_path(data = streams_gg, 
            aes(x = long, y = lat, group = group), 
            color = "lightblue", size = 1.5)
ggsave("shp_and_ggmap.png", dpi = 72) 

Static maps - with data and formatting

library(dplyr)
alum_data <- plot_data %>% filter(Analyte == "ALUMINUM")
g + 
  geom_path(data = streams_gg, 
            aes(x = long, y = lat, group = group), 
            color = "lightblue", size = 1.5) +
  geom_point(data=alum_data, 
             aes(x = Longitude, y = Latitude, size = Result, fill = Result),
             shape = 21, alpha = 0.9) +
  scale_fill_gradient(low = "#1f77b4", high = "#d62728") +
  scale_size(range = c(2, 8)) +
  guides(size = FALSE) +
  coord_equal() 
ggsave("static_map.png", dpi = 72) 
head(alum_data)
##    Analyte    Result Units Latitude Longitude
## 1 ALUMINUM  9120.859 mg/kg 36.59625 -109.1736
## 2 ALUMINUM  7429.408 mg/kg 36.59625 -109.1736
## 3 ALUMINUM  1657.721 mg/kg 36.61342 -109.1356
## 4 ALUMINUM  1709.072 mg/kg 36.61342 -109.1356
## 5 ALUMINUM  5867.922 mg/kg 36.56392 -109.2120
## 6 ALUMINUM 46050.939 mg/kg 36.56392 -109.2120

Static map - with data and formatting

Using R in place of ArcGIS

3D maps

Lollipop example, using scatterplot3d

library(scatterplot3d)
with(alum_data, {
  lollipop_plot <- scatterplot3d(streams_gg$long, streams_gg$lat, 
                                 rep(0, length(streams_gg$long)),  # x y and z axis
                              color = "blue", pch = 16, type = "p",
                              angle = 120,            
                              scale.y = 1.75,           
                              main = "Example Lollipop plot",
                              xlab = "Longitude",
                              ylab = "Latitude",
                              zlab = "Concentration values",
                              xlim = c(-109.273, -109.15),
                              ylim = c(36.50, 36.59),
                              zlim = range(0, 46000),
                              grid = TRUE,
                              box = FALSE)
  
  # add the legend
  legend("topright", inset = 0.07,      
         bty = "n", 
         title = "Type of Results",
         c("High","Low"), 
         fill = c("#1B9E77", "#D95F02"))
  
  # add the lollipop points
  lollipop_plot$points3d(Longitude, Latitude, Result,        
                      col = "#282830", 
                      pch = 21, 
                      bg = ifelse(Result > 10000, "#1B9E77", "#D95F02"), 
                      lwd = 1,        
                      type = "h",      
                      cex = (Result / 50000) + 1)
})

Lollipop example

Lollipop example

lollipop plot

lollipop plot

(Generated using a modified version of the scatterplot3d function: https://github.com/USEPA/R-User-Group/tree/master/contributedCode/scatterplot3dMap)

Reference: Beaulieu, et al. 2016. Estimates of reservoir methane emissions based on a spatially balanced probabilistic-survey. Limnology and Oceanography, 61: S27-S40.

3D map example (interactive)

library(rgl)
library(dplyr)
plot_3d <- with(alum_data, 
  plot3d(Longitude, Latitude, log(Result),        # x y and z axis
         col = ifelse(Result > 10000, "#1B9E77", "#D95F02"), size = 5,
         type = "p"))
rglwidget(elementId = "plot3drgl") # to show in presentation

3D map- plotly

## plotly recommends the development version of ggplot2, but will also work with the
## latest version of the package, version 2.2.1 
# library(devtools)
devtools::install_github('hadley/ggplot2')
library(ggplot2)
library(plotly)
p <- plot_ly(alum_data, 
             x = ~Longitude, y = ~Latitude, z = ~Result,
             marker = list(color = ~Result, 
                      colorscale = c('#FFE1A1', '#683531'), showscale = TRUE)) %>%
  add_markers() 

3D map- plotly

p

Interactive maps- plotly

interactive_map <- g +  
  geom_point(data = alum_data, 
             aes(x = Longitude, y = Latitude, size = Result, fill = Result),
             color = "#000000", shape = 21, alpha = 0.8) +
  scale_fill_gradient(low = "#1f77b4", high = "#d62728",
                      guide = "legend") +
  scale_size(range = c(2, 8)) +
  coord_equal() +
  ggtitle("ALUMINUM")

Interactive maps- plotly

ggplotly(interactive_map + theme_void(), filename="plotly")

Interactive maps cont’d

Animated maps

library(animation)
ani.record(reset = TRUE)
for(a in unique(plot_data$Analyte)) { 
  plot_data_a <- subset(plot_data, Analyte == a)
  animated_maps <- g + 
    geom_path(data = streams_gg, 
              aes(x = long, y = lat, group = group),
              size = 1.5, color = 'lightblue') +
    geom_point(data = plot_data_a, 
               aes(x = Longitude, y=Latitude, size = Result, fill = Result),
               shape = 21, alpha = 0.8) +
    scale_fill_gradient(low = "#1f77b4", high = "#d62728") +
    scale_size(range = c(2, 8)) +
    guides(size = FALSE) +
    coord_equal() +
    ggtitle(a)
  
  print(animated_maps)
  ani.record()
}
oopts = ani.options(interval = 1)
ani.replay()
saveHTML(ani.replay(), img.name = "animation_plot", 
         htmlfile = "animation_5metals.html", ani.width=600, ani.height=400)

Animated maps

Spatial data exploration- Moran’s I

# Show neighbors on map of subbasins
library(sp)
class(lulc.sp1)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
spplot(lulc.sp1, zcol = "PCTHERB", col = "black", 
       main = "Percent Herb")

Spatial data exploration- Moran’s I

library(spdep)

## determine who shares a border
ctchmt.nb1 <- poly2nb(lulc.sp1, queen = TRUE)
# put neighbors into list
ctchmt.nbwts.list <- nb2listw(ctchmt.nb1, style = "B")

moran.test(lulc.sp1$PCTHERB, ctchmt.nbwts.list)
## 
##  Moran I test under randomisation
## 
## data:  lulc.sp1$PCTHERB  
## weights: ctchmt.nbwts.list  
## 
## Moran I statistic standard deviate = 2.9826, p-value = 0.001429
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.212182899      -0.016949153       0.005901866

Spatial data exploration- Moran’s I cont’d

par(mar=c(5, 5, 2, 2))
moran.plot(lulc.sp1$PCTHERB, listw = ctchmt.nbwts.list, 
           labels = lulc.df1$PCTHERB,
           xlab = "PCTHERB", ylab = "Spatially Lagged PCTHERB",
           main = "Moran Plot for PCT HERB")

Additional Packages

Appendix: Interactive network graph

## The network graph may require the following graph package that need to be 
## installed from source.  (This package is not generally required to create  
## similar graphs; they are only required for the example.  See the tutorial for  
## other examples based on csv files.)
# source("https://bioconductor.org/biocLite.R")
# biocLite("graph")
# biocLite("Rgraphviz")

library(magrittr)  ## for %<>% function
library(visNetwork)  ## to create the interactive network graph
library(jsonlite)
library(igraph)
network_interactive <- function(g, nodes, coords){
  ## first, switch graph construction to data frames:
  nNodes <- length(nodes)
  nodes_df <- matrix(NA, nrow=nNodes, 4)  # save first 4 attributes
  for (i in 1:nNodes){
    nodes_df[i, ] <- unlist(jsonlite::fromJSON(nodes[[i]]))[1:4]
  }
  nodes_df <- as.data.frame(nodes_df, stringsAsFactors=FALSE)
  
  ## because edges uses the names, need id to equal names of nodes, not numbers
  names(nodes_df)[1:3] <- c("old_id", "id", "type")
  
  ## and id column needs to go first: (only for igraph package, not visNetwork)
  ## igraph packge was used to help decide node coordinates
  nodes_df <- nodes_df[, c('id', 'old_id', 'type')]

  g_edges <- Rgraphviz::buildEdgeList(g)
  nEdges <- length(g_edges)
  edges_df <- NULL
  for (i in 1:nEdges){
    tmp <- c(from=g_edges[[i]]@from,
             to=g_edges[[i]]@to,
             unlist(g_edges[[i]]@attrs))
    edges_df %<>% bind_rows(tmp)
  }
  edges_df$weight <- as.numeric(edges_df$weight)
  
  ## Add layer info/ color column
  nodes_df$col_layer <- 
    ifelse(nodes_df$id %in%  grep("MainInput", nodes_df$id, value=T), 1, # 
           ifelse(nodes_df$id %in% grep("OtherInput", nodes_df$id, value=T), 3,
           ifelse(nodes_df$id %in%  grep("Midvalue", nodes_df$id, value=T),  4,
           ifelse(nodes_df$id %in%  grep("MidEqn", nodes_df$id, value=T),  5,
                  6))))
  
  ##  add coordinates to determine layout of each node:
  nodes_df <- full_join(nodes_df, coords)

  ## specify colors to use for each col_layer
  graph_colors <- c("gold", "darkorange", "tomato", 
                    "palegreen", "seagreen", "royalblue")
  # frame = borders
  frame_colors <- c("gold3", "darkorange3", "tomato4", 
                    "palegreen3", "seagreen4", "royalblue4")
  ## (visNetwork doesn't like color names with numbers)
  frame_colors_rgb <- rgb(t(col2rgb(frame_colors)), maxColorValue = 255)

  ## Add attributes so that the graph looks good:
  nodes_df$shape <- "ellipse"
  nodes_df$color.background <- graph_colors[nodes_df$col_layer] 
  nodes_df$color.border <- frame_colors_rgb[nodes_df$col_layer]
  nodes_df$color.highlight.border <- "darkred"
  
  edges_df$arrows <- "to"  # draw arrows on the edges.
  edges_df$color.highlight <- "black"
  
  network <- visNetwork::visNetwork(nodes_df, edges_df, width="100%", physics = TRUE) %>%
    visNetwork::visIgraphLayout() %>%
    visNetwork::visOptions(highlightNearest = TRUE, nodesIdSelection = TRUE,
               manipulation = FALSE) %>%
    visNetwork::visEdges(color = "slategray", smooth = list(roundess = 1)) %>%
    visNetwork::visNodes(font = list(size = 22), shape="ellipse")

  return(network)
}

Appendix: Interactive network graph

network_interactive(g = networkDag, nodes = networkNodes, coords = networkCoords)

Appendix: Color

Appendix: Color

par(mar=rep(0, 4))  ## change margins so plot fills entire space
## default colors:
plot(1:8, pch=16, cex=5, col=1:8)

## Defining a new palette
new_palette <- c("darkred", "chartreuse", "turquoise", "purple",  
                 "gray45", "plum", "black", "#F08800")
palette(new_palette)
plot(1:8, pch=16, cex=5, col=1:8)

palette("default")  ## return palette to default colors

Version Control

sessionInfo()
## R version 3.4.2 (2017-09-28)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Sierra 10.12.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] igraph_1.1.2       jsonlite_1.5       visNetwork_2.0.1  
##  [4] magrittr_1.5       spdep_0.6-15       Matrix_1.2-11     
##  [7] bindrcpp_0.2       plotly_4.7.1       dplyr_0.7.4       
## [10] rgl_0.98.1         ggplot2_2.2.1.9000 rgdal_1.2-13      
## [13] sp_1.2-5          
## 
## loaded via a namespace (and not attached):
##  [1] httr_1.3.1          maps_3.2.0          tidyr_0.7.1        
##  [4] viridisLite_0.2.0   splines_3.4.2       gtools_3.5.0       
##  [7] shiny_1.0.5         assertthat_0.2.0    expm_0.999-2       
## [10] stats4_3.4.2        yaml_2.1.14         LearnBayes_2.15    
## [13] backports_1.1.1     lattice_0.20-35     glue_1.1.1         
## [16] digest_0.6.12       colorspace_1.3-2    htmltools_0.3.6    
## [19] httpuv_1.3.5        plyr_1.8.4          pkgconfig_2.0.1    
## [22] devtools_1.13.3     gmodels_2.16.2      purrr_0.2.3        
## [25] xtable_1.8-2        scales_0.5.0.9000   gdata_2.18.0       
## [28] jpeg_0.1-8          ggmap_2.6.1         git2r_0.19.0       
## [31] tibble_1.3.4        withr_2.0.0         BiocGenerics_0.22.1
## [34] lazyeval_0.2.0      proto_1.0.0         mime_0.5           
## [37] deldir_0.1-14       memoise_1.1.0       evaluate_0.10.1    
## [40] nlme_3.1-131        MASS_7.3-47         graph_1.54.0       
## [43] tools_3.4.2         data.table_1.10.4-1 geosphere_1.5-5    
## [46] RgoogleMaps_1.4.1   stringr_1.2.0       munsell_0.4.3      
## [49] compiler_3.4.2      rlang_0.1.2         grid_3.4.2         
## [52] rjson_0.2.15        htmlwidgets_0.9     crosstalk_1.0.0    
## [55] base64enc_0.1-3     labeling_0.3        rmarkdown_1.6      
## [58] boot_1.3-20         gtable_0.2.0        curl_3.0           
## [61] reshape2_1.4.2      R6_2.2.2            knitr_1.17         
## [64] bindr_0.1           rprojroot_1.2       Rgraphviz_2.20.0   
## [67] stringi_1.1.5       parallel_3.4.2      Rcpp_0.12.13       
## [70] mapproj_1.2-5       png_0.1-7           coda_0.19-1